Classical dimers with aligning interactions on the square lattice 
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We present a detailed study of a model of close-packed dimers on the square lattice with an 
interaction between nearest-neighbor dimers. The interaction favors parallel alignment of dimers, 
resulting in a low-temperature crystalline phase. With large-scale Monte Carlo and Transfer Matrix 
calculations, we show that the crystal melts through a Kosterlitz-Thouless phase transition to give 
rise to a high-temperature critical phase, with algebraic decays of correlations functions with expo- 
nents that vary continuously with the temperature. We give a theoretical interpretation of these 
results by mapping the model to a Coulomb gas, whose coupling constant and associated exponents 
are calculated numerically with high precision. Introducing monomers is a marginal perturbation at 
the Kosterlitz-Thouless transition and gives rise to another critical line. We study this line numer- 
ically, showing that it is in the Ashkin- Teller universality class, and terminates in a tricritical point 
at finite temperature and monomer fugacity. In the course of this work, we also derive analytic 
results relevant to the non-interacting case of dimer coverings, including a Bethe Ansatz (at the free 
fermion point) analysis, a detailed discussion of the effective height model and a free field analysis 
of height fluctuations. 

PACS numbers: 05.20.-y, 05.50.+q, 64.60.Cn, 64.60.Fr, 64.60.Kw 



I. INTRODUCTION 

The problem of lattice coverings by "hard" objects, 
and dimers in particular, is ubiquitous in classical statis- 
tical mechanics. The formulation of the problem of dimer 
coverings goes back to the thirties 0. The combinato- 
rial problem of finding the exact number of such coverings 
has been solved in the early sixties for planar 2d lattices 
by means of Pfaffian techniques 0, , which have been 
extended to calculate dimer-dimer and monomer (i.e. un- 
paired sites)-monomer correlation functions ,4j . Notwith- 
standing the intrinsic mathematical beauty of this prob- 
lem it also plays a central role in statistical physics 
due to its relationship to Ising or height models 0. 
Dimer coverings of bipartite graphs in three dimensions 
have also been recently shown to be connected to gauge 
theories 0. Dimer models have also recently regained 
interest because Quantum Dimer Models (QDM), origi- 
nally introduced by Rokhsar and Kivelson [8j , are among 
the simplest systems which exhibit ground-states with 
topological order and fractionalization 0. 

In this work, we study a model of interacting clas- 
sical dimers on the square lattice, with an interaction 
that favors dimer alignment. The dimer coverings are 
close-packed, i.e. there are no sites left uncovered by a 
dimer (monomers). We now describe the plan of the pa- 
per: we first introduce the model and its simple limits in 
Sec. [n] In Sec. IIIII we introduce the transfer matrix of 
the model and describe how its critical exponents in the 
non-interacting (infinite temperature) limit can be red- 



crived by the Bethe Ansatz technique. Unfortunately, 
the interacting model does not seem to be integrable by 
a straightforward extension of this approach. We there- 
fore go on, in Sec. IIVI to describe two complementary 
numerical simulation schemes: Monte Carlo (MC) and 
Transfer Matrix (TM) calculations. The results of the 
MC simulations are presented in Sees. [V] IVII and IVIII 
we find that the model possesses a low-temperature crys- 
talline phase separated by a Kosterlitz-Thouless (KT) 
transition |lfjj from a high-temperature critical phase 
with floating exponents. We account for all these findings 



in Sec. IVIIII where we give a theoretical interpretation 
in terms of a Coulomb gas (CG) picture ^lj. This map- 
ping moreover allows use to make specific predictions on 
the high-temperature phase that are successfully tested 
with high-precision TM and MC calculations. The CG 
description implies that the introduction of monomers is 
a marginal perturbation at the KT point and hence leads 
to the emergence of another critical line. We study this 
numerically and find it to be in the Ashkin- Teller uni- 
versality class; the line terminates in a tricritical point 
at finite temperature and monomer fugacity. The finally 
obtained phase diagram is presented in Fig.Q] We finally 
discuss the connections to other models in classical sta- 
tistical physics and the implications of our findings for 
quantum models in Sec. IIXI and conclude in Sec. A 
short account of the results presented here was given in 
Ref.[H 
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FIG. 1: (Color online) Phase diagram of the interacting dimer 
model in the temperature T, monomer fugacity £ plane (see 
text for definitions). The solid lines represent second-order 
phase transition lines with continuously varying exponents. 
When no monomers are allowed (£ = 0), the first critical line 
terminates at T c = 0.65(1) and separates the high-T criti- 
cal phase from a long-range order crystalline phase through a 
Kosterlitz-Thouless phase transition. Allowing for monomers 
(£ 0) creates the second critical line separating the low 
T crystalline phase from a monomer-dimer (massive) liquid 
phase at high T. This line terminates in a multicritical point 
at T* = 0.29(2), where it changes nature to become a first 
order line (dashed line). Simple energetic arguments (see 
Ref. H2T) predict that the first order transition temperature 
scales as 1/(2 ln(£)) when £ — » 00. 



II. THE MODEL 

We study a model of interacting close-packed dimcrs 
on the square lattice, defined in the following way 



c 



(1) 



The sum in the partition function Z is over all fully- 
packed dimer coverings of the square lattice c. To each 
dimer covering c, we assign the energy E c which sim- 
ply counts the number 7V c (z) + N c (l l) of plaquettes with 
parallel (horizontal or vertical) dimers in the covering c. 
\v\ = 1 sets the energy scale (T is the temperature). The 
sign of v determines the nature of the interactions be- 
tween the nearest-neighbor dimers : v < correspond to 
aligning interactions between dimers, v > favors con- 
figurations with staggered occupation of dimcrs. In this 
paper, we will consider v = — 1, the so-called columnar 
case. 



The model is illustrated in Fig. [21 where a dimer cover- 
ing of the lattice is represented, and the plaquettes con- 
tributing a factor +v to the energy of this configuration 



FIG. 2: Illustration of the interacting dimer model : we con- 
sider dimer coverings of the square lattice where each pla- 
quette (marked with a cross) with a pair of parallel nearest 
neighbor dimers contributes +v to the energy (the energy of 
this dimer covering is 7v). 



FIG. 3: The four columnar ground states. 



are identified by a cross. It is straightforward to see that 
at zero temperature T = 0, the configurations that mini- 
mize the energy are the four states represented in Fig. |3 
where the dimers are aligned in columns. This four-fold 
degenerate ground state spontaneously breaks transla- 
tion and 7r/2-rotational symmetries. The first excitation 
above these ground-states are obtained by flipping two 
parallel dimers around a plaquette; the system has a gap 
(it costs a finite energy 2v to flip the two dimers) and 
the columnar order is therefore expected to subsist at 
(possibly small but) finite temperature. 



On the other hand, the partition function at infinite 
temperature T — 00 is simply the unweighted sum over 
all possible dimer coverings of the square lattice,and the 
model can be solved exactly at this point 0, The 
T = 00 point is critical, with correlation functions dis- 
playing an algebraic dependence with distance Q : dimer- 
dimer correlation functions decay as 1 /r 2 and monomer- 
monomer correlation functions as 1/ ^/r for large distance 
r. We postpone the precise definitions of these correla- 
tion functions to Sec. IIIII below where we rederive the 
results for the critical exponents using another exact ap- 
proach, the coordinate Bethe Ansatz. Although we have 
not been able to solve the interacting dimer problem (fi- 
nite temperature), the Bethe Ansatz technique can po- 
tentially go beyond free fermion problems (contrary to 
the Pfaffian methods of Refs.HHH. 

The Bethe Ansatz method also serves to illustrate that 
the critical nature of the dimer covering problem is inti- 
mately linked to the bipartite nature of the square lattice 
(non-bipartite lattices present a dimer liquid behavior 
with a finite correlation length 0). Unfortunately, the 
introduction of interactions appears to break the integra- 
bility of the model. 
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We end up the introduction with some historical notes 
on this model. The model Eq. Q was first introduced 
in the physics of liquid crystals [ljj and not developed 
further in this context to our best knowledge. This is 
likely due to the fact that the quest was there to look 
for microscopic models where a true liquid crystal phase 
exists, and not a crystalline state such as the one dep icted 
in Fig. [3J Later on, Brankov and coworkers [1J| also 
studied the same model with Monte Carlo methods but 
missed the true critical behavior of this problem. In a 
recent publication [l^l we described the physics of the 
undoped model and sketched the existence of a critical 
line with central charge c = 1 ending at a tri-critical point 
at finite doping. This has then been followed by other 
studies, including construction of quantum models [T^L 
IT^ I with ground-state wave- functions described by the 
partition function in Eg • JM ; further investigations of 
the doped monomer case Il7| and generalization 

to three-dimensional lattices llSl. 



III. NON-INTERACTING DIMERS AS FREE 
FERMIONS 



Pi 



m 



FIG. 4: The row-to-row transfer matrix 



2 3 4 5 6 1 



1 2 3 4 5 6 

(a) 

2 3 4 5 6 1 



2 3 4 

(c) 



2 3 4 5 6 1 



1 2 3 4 5 6 

(b) 

2 3 4 5 6 1 



1 2 3 4 5 6 
(d) 



We study dimer coverings (each site is paired with ex- 
actly one of its neighbors) of the square lattice. In this 
section (and only in this section), we give a fugacity tu 
to each horizontal dimer, and 1 to vertical ones. This 
is the model of non-interacting dimers, solved by com- 
binatorial methods 0, 0, 3 ■ We will introduce the TM 
of this model and we show how to compute the partition 
sum and the correlation functions by the Bethe Ansatz 
method. 



FIG. 5: Rules for the dynamics of the particles on a strip of 
width L — 6. Particles are represented by zigzag lines. The 
rules are different for particles starting from an even site (a, 
b, c) and an odd site (d). 



B. Conservation law 



A. Transfer matrix 

The partition sum of the model is: 

^7 ^ ^ ^^£horizontal dimers 

dimer config. 

On a strip of width L, we define the state of a row as the 
"occupation numbers" a = (oti, . . . , «l) of the vertical 
edges, where on is equal to 1 if the i-th vertical edge is 
occupied by a dimer, and otherwise. There are 2 L con- 
figurations of a row, so Z can be written as the trace of 
a 2 L -dimensional transfer matrix T. Wc impose periodic 
boundary conditions (PBC), i.e. the index i is consid- 
ered modulo L. Given two line configurations a and /?, 
the matrix element Tp a is the sum of the Boltzmann 
weights associated with the horizontal dimer configura- 
tions [i compatible with a and j3: 

T Pa = Yl (3) 

fl(«,/3) 

as illustrated in Fig. 0] The compatibility criterion 
n\(a,P) can be expressed formally as follows: 

Vie {!,..., L} : m + + cti + /3 i+1 = 1. (4) 



For convenience, we introduce a shift at each row in 
the numbering of columns (see Fig. Wc call particle 
an empty even vertical edge or an occupied odd vertical 
edge. Let us show that the number of particles is con- 
served, and let us give at the same time the rules for the 
dynamics of the particles (see the corresponding Fig. [5j|. 

Particle on an even vertical edge: If a particle sits 
on the vertical edge a 2 j, then a 2 j = 0. The site above 
this edge must be visited once, so one of the variables 
M2ji 02j+i, M2j+i must be equal to 1. In the first (resp. 
third) case, this implies that (3 2 j (resp. fi 2 j+2) is zero. 
Therefore, in the next row, there is a particle on the 
edge faj, foj+i or /3 2 j+2- 

Particle on an odd vertical edge: If a particle sits on 
the vertical edge OL2j—u then Oi2j-i = 1- The site above 
this edge has already been visited, so @2j must be zero: 
in the next row, the particle sits on the edge foj- 

The TM is block-diagonal, each block representing a 
sector with fixed number of particles n. We call T'™- 1 the 
TM block in the n-particle sector. Note that the lattice 
width L must be even, because, for an odd lattice width 
with PBC, the number of particles is not conserved. 
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C. One-particle sector 

The action of T on a one-particle state $ is: 

(T$) (2j) = w$(2i-2) + $(2j-l)+cj$(2j) 
(T$) (2j + l) = 4>(2j) 

. (5) 

We want to take advantage of the translational invariance 
to diagonalize T^K Define the two-step cyclic permuta- 
tion J of the sites by its action on a one-particle state 



(J$)(x) = <f>(x 



(6) 



On a lattice of even width with PBC, the operator J 
commutes with If z satisfies the condition z L = 1, 

the eigenspace of J with eigenvalue z 2 is generated by 
the two vectors 



**(2i) 

**(2j) 



2j 



*,(2i-l) = 
$ 2 (2j-l) = z 2 ^ 1 



(7) 



Note that $_ 2 = $ 2 and I>_ 2 = -¥ z . 

More generally, let z be a complex number of modulus 
unity. The block of in the basis ($ 2 , $ 2 ) is: 



w (1 + z~ 2 ) z" 







(8) 



This matrix has eigenvectors V'z , ^ with the eigenvalues 
A(z),A'(z) satisfying : 



A(z) + A'(z) = w (1 + z" 2 ) 
A(z) A'(z) = -z- 2 

One can then write: 

z = cxp (ik) 
A(z) = exp[h + i((t> + 9)] 
A'(z) = exp[-h + i((j>-6)] 



(9) 



(10) 

(11) 
(12) 

A(z) 



where k, 4>, 9 are real and h is nonnegative. 
is the eigenvalue with greatest modulus. With this 
paramctrization, Eqs. © imply: 



cosh(2/i) = 1 + 2uj 2 cos 2 k 



(13) 



when |z| = 1. Recalling that cosfc > 0, this relation can 
be inverted and we obtain: 



h(k) = log u> cos k + (1 + u! 2 cos 2 k)~- 



(14) 



when — 7r/2 < k < tt/2. A useful quantity for the com- 
putation of finite-size effects is h'(w/2) = —uj. 



D. Two-particle sector, scattering amplitude 



Consider the action of the TM on the two-particle vec- 



tor: 



tpnixi, x 2 ) = ip zi (x 1 )ip Z2 {x 2 ) , xi < x 2 (15) 



The TM changes the positions of the particles from 
(xi,x 2 ) to (2/1,2/2)- For fixed positions 2/1 < y 2 , let us 
look at the initial states leading to these positions: 

• if 2/2 > 2/1 + 2 or (2/1,2/2) = (2j - 1, 2j + 1), then for 
each particle all initial states arc allowed, thus: 

(Tipi 2 )(yi,y 2 ) = ^2 ^zi(x 1 )i) z . 2 (x 2 )T yi . Xl Ty^ 

X1,X 2 

= A(zi)A(z 2 )il) Zl (yi)^ Z2 {y 2 ) 

• if (2/1,2/2) = (2j — 1, 2j), all initial states are allowed 
except x\ = x 2 = 2j — 2. One has to subtract the 
corresponding term in the action of the matrix T: 

(T^ia)(2j-l,2j) = A(zi)A(z 2 )^ 1 (2. 7 -l)^ 2 (2j) 
-u>i) zl (2j-2)i> Z2 {2j-2) 

• if (2/1,2/2) = (2i,2j + 1) or (2/1,2/2) = (2j,2j + 2), 
all initial states are allowed except x\ = x 2 = 2j. 
Similarly to the previous case: 

(2Via)(2j,2j + l) = A(zi)A(z 2 )^ 1 (2j>. 2 (2. 7 + 1) 
-LJiP zl (2j)^ Z2 (2j) 

(2> 12 )(2j,2j + 2) = A(z 1 )A(z 2 )V Zl (2j>. 2 (2. 7 + 2) 
-LJtP zl (2j)^ Z2 (2j) 

Note that all the interaction terms in T^i2 are symmetric 
functions of the momenta ki,k 2 . As a consequence, the 
antisymmetric combination: 

ip(xi,x 2 ) = ip Zl {xi)'ip Z2 (x 2 ) - ip Z2 (xi)ip zl (x 2 ) (16) 

is an eigenvector of the matrix T with eigenvalue 
A(zi)A(z 2 ). 



E. Periodic boundary conditions, position of the 
solutions 



The analogous construction in the n-particles sector 
gives the eigenvectors and eigenvalues: 

ip(x 1: ...,x n ) = ^2 e (- p )V , 2 P1 (xi)... tp Zpn (x n ) (17) 



A{z l ,...z n ) = A( Zl ) . . . A(z n ) 



(18) 



where the sum is over all permutations of the integers 
1, . . . , n and e(P) is the signature of the permutation P. 
PBC yield: 

Vj = (-I)"" 1 (19) 

The solutions of these equations lie on the unit circle, 
which justifies the discussion in section IIII CI The mo- 
menta kj are given by: 



k o = T" 7 ^ J J" e Z n odd 

k j = T7 & + I ) . h 6 z n even 



(20) 
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n odd 



n even 



FIG. 6: Position of the vacancies for a strip of width L = 12. 
Vacancies are represented by black points on the unit circle, 
in the complex z plane. 



See Fig. for a graphical representation of the vacan- 
cies on the unit circle. For any z on the unit circle, the 
TM has an eigenvalue A with modulus greater than one. 
Therefore, the number of particles that maximizes the 
total eigenvalue of T is either the greatest even value or 
the greatest odd value for n. Since the fcj's are distinct, 
lie in the interval [— 7r/2,7r/2] and are spaced by 2ir/L, 
the maximum number of particles is L/2. 



F. Thermodynamic limit 

In this section, the discussion is restricted for simplicity 
to a system of width multiple of four: L — Ap. According 
to the previous section, the leading sector is defined by 
the greatest eigenvalue either in the sector n = 2p or in 
the sector n = 2p — 1. As will be shown in a few lines, 
the correct choice for the leading sector is n = 2p. 

For a system of finite width and infinite length, the free 
energy density per surface unit in the n-particles sector 
is defined by: 



An) 
■I L 



-1okA (?i) 

tug i v ma; 



L 



(21) 



where A^ax is the eigenvalue of T with greatest modulus 
in the n-particles sector. The corresponding quantity for 
L = Ap and n = 2p is: 



,(L/2) 
J L 



2 p_1 



(22) 



where the function h(k) is given by Eq. 1|14|) . When L 
goes to infinity, this quantity tends to the limit /oo: 

1 r r 

foo = — fog lu cos k+ (l +w 2 cos 2 fc)2 dk (23) 

71" JO 1 J 

in agreement with formula (17) of Ref. 2} In the isotropic 
case lu = 1: 



/oo(w = l) = 



2G 



where G is the Catalan constant: 



G 



(24) 
(25) 



The asymptotic behavior of /| L ^ 2 ' 
Euler-Maclaurin formula: 



is derived from the 



JL/2) 
J L 



6L 2 



(26) 



We expect the critical point to have conformal symmetry 
in the isotropic case, with a central charge c = 1. If one 
particle is removed {n = 2p — 1), the solutions Zj all get 
shifted (see Fig.|HJ). The Euler-Maclaurin formula yields: 



with 



/r-^/oo-^'^-^^+o^ 2 ) (27) 



1(e) = ' h(k) dk = -ti ^e 2 /2 + o(e 2 ). (28) 



Finally we obtain: 

AL/2-1) _ AL/2) 
JL — JL 



2L^ 



(29) 



This proves that in the thermodynamic limit the leading 
sector is indeed given by n = 2p. In the isotropic case, 
the critical exponent corresponding to the removal of one 
particle is A"i = 1/4 (see definition and discussion in 
Sec lIVBl below^l. Now if two particles are removed from 
the leading sector, the other zfs are not shifted, because 
the number of particles remains even. The only effect is 
a decrease of free energy caused by the absence of the 
two particles: 



AL/2-2) 
JL 



/i i/2) -|- + (L-) 



(30) 
(31) 



In the isotropic case, the critical exponent corresponding 
to this process (see again Sec. lTVRl helow'l is X 2 = 1. 



G. Introducing interactions 

A TM for the interacting model Eq. 0} can be written 
down by generalizing the working of Sec. IIII Al To this 
end, the basis states (a) must encode not only the oc- 
cupation numbers of a row of vertical edges (as before), 
but also the occupation numbers of the preceding row of 
horizontal edges, as shown in Fig. 

The transfer matrix T is most easily defined by giving 
its sparse matrix decomposition 



L-l 



T 



T 2 Ti, T k =Y[T, 



(i) 
fc ' 



(32) 



where the matrix encodes the interactions at pla- 

(i) 

quette i and T 2 imposes the dimer constraint at vertex 
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FIG. 7: The row-to-row transfer matrix for the interacting 
case. 



i. More precisely, evolves a 2 i+i into (3 2i+ i and has 
matrix elements (we set W = e~ v / T ) 

T^\a 2i a 2 i + i\p 2 i + ia 2i+2 ) = W{a 2i a 2 i +2 + a^i+ife+i) , 

(i) 

whereas T 2 evolves a 2 i into p 2 % and has matrix elements 

We have attempted to diagonalize T using the Bcthc 
Ansatz method, using a straightforward generalization of 
the working exposed in the preceding subsections but we 
failed to obtain a consistent determination of the scat- 
tering amplitudes S(zi,Zj). Most likely, this means that 
the interacting dimcr model is not integrable. In the re- 
mainder of the paper, we therefore study the model using 
numerical and (non-rigorous) field theoretical methods. 



IV. NUMERICAL METHODS 

Our numerical methods consist in MC simulations and 
exact diagonalization of the TM. We now describe these 
two methods in turn. 

A. Monte Carlo calculations 

We use a MC directed- loop (or directed- "worm" ) algo- 
rithm 0| . This method allows to make non-local moves 
in the dimer configurations by changing the positions of 
dimers along a closed loop, which can be quite large. This 
results in small autocorrelation times in the MC process, 
and permits to treat large systems (up to 512 x 512 in 
this study). Moreover, the directed-loop algorithm cap- 
tures the physics of test defects (monomers) in the dimcr 
configuration as we discuss below. The algorithm indeed 
allows to calculate monomer-monomer correlation func- 
tion; conversely, this indicates that the performance of 
the algorithm is dictated by the physical properties of 
test monomers in the different physical phases. 

For sake of completeness, we briefly describe below the 
algorithm following Rcf. [l|| and specifying minor details 
where our specific implementation differs. One MC sweep 
of the algorithm consists of the three following steps: 



1. The worm, which can be seen as constituted by two 
monomers (head and tail), is initially placed on top 
of a dimer configuration at a random site i = iq. 

2. The site i is connected to a neighboring site j by a 
dimer in the background configuration (this dimer 
is noted The head of the worm is moved to 
j and the dimer (i,j) is removed, leaving the site 
j with no dimer attached to it. Out of the four 
neighbors of j, one (which we call k) is selected ac- 
cording to a local detailed balance rule (see below). 
A dimer is put between j and k. 

3. If k = io, the worm is finished and we are left with 
a new valid dimer configuration. Otherwise, we 
rename i — k and go back to step 2. 

How does the worm, sitting at site j (and coming from 
site i), choose the site k where a dimer will be put in 
step 2 ? For this, we consider the weights (respec- 
tively W(jk)) contributed to the partition function by a 
dimer located between sites i and j (respectively j and 
k). In the model of Ref. H each dimer is given a cer- 
tain fugacity and thus contributes solely a certain weight 
to the partition function. In our model, the weight of 
a dimer is given by ttffe-) = exp(—v.Nij/T) where 

Nij € {0, 1, 2} is the number of nearest neighbors parallel 
to the dimer Once these weights are known, the 

probability P((i,j) — > (j,k)) to select a given site k is 
imposed to satisfy a local detailed balance rule: 

P((iJ) (i>*0)w(ij) = P(U,k) -> (i,j))w {jk) . (33) 

This leads to a set of equations ( "directed-loop" equa- 
tions l|j) corresponding to all the possibles values of 
local dimer configurations and the corresponding num- 
bers Nij. These equations are underdetermined, and we 
impose by experience p(l l2l| to minimize the bounce 
processes P((i,j) — ► i.e. the case where the site 

k is chosen to be the origin site i (the worm backtracks 
in its own path, which is a priori quite useless). For 
the specific model of Ref. fill a solution minimizing the 
bounce probabilities and satisfying the local detailed bal- 
ance equation was found analytically. More generally, 
such a solution can always be found numerically with 
linear programming techniques [2l| ■ Please note that at 
T = oo, the worm simply performs a random walk in the 
dimer configuration (more precisely, all the even steps in 
the walk are purely random, the odd ones are dictated 
by the underlying dimer configuration). 

Taking a snapshot of the configuration during the 
worm construction shows that two test monomers have 
been inserted in the dimer configuration, and thus con- 
nect the behavior of the worm to the monomer corre- 
lation function. The fact that the worm walk is locally 
detailed balance actually imposes the histogram of the 
distance r between the worm's head and tail to be propor- 
tional (up to a small correction factor) to the monomer- 
monomer correlation function M (r) (see precise defini- 
tion in Sec. IVIIf) . The proof of this statement can be 
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worked out along the lines of Ref. OH]. The only subtlety 
is the following: the measurement of the correlation func- 
tion is made at step 2, before the selection of the next 
site k. Since all the future dimer positions (j, k) are not 
equivalent (they will contribute differently to the parti- 
tion function) and since the next dimer position is not 
yet decided, we have to correct the monomer-monomer 
correlation estimator by the inverse of the total weight 
contributed by all possible future positions, i.e. we in- 
crement the estimator of M(r) (with r the position dif- 
ference vector between sites i and j) by Wf 1 where 
Wj = J2k w U k )- ^ a ^ future position dimers are equiva- 
lent (as in the model of Ref. fl9f) . this factor is constant, 
and we can just simply identify the histogram of the dis- 
tance r between the worm's head and tail to M(r). 

The worm algorithm therefore possesses the nice fea- 
ture of being able to calculate M(r), even if monomers 
are not allowed in the model. This also indicates that 
the worm algorithm performances is bound to follow the 
physics of monomers: if the monomers are confined, the 
worms will be short, resulting in a merely local algorithm 
- which is known to display poor performances (for ex- 
ample ergodicity problems). If the monomers are decon- 
fined, worms will be long and will update a massive num- 
ber of dimers - resulting in small autocorrelation times. 

The technical details of the MC calculations are as fol- 
lows: simulations were performed on N — L x L samples, 
up to L = 160 for the full T range, and up to L = 512 
for correlation functions for a few chosen temperatures. 
PBC are assumed. Each MC sweep is constituted by a 
number of worms such that all the links of the lattice are 
visited once on average by a worm. For each parameter 
set, a total between 10 6 and 10 7 sweeps was performed. 

B. Transfer matrix calculations 

The TM for the interacting dimer model was defined 
above in Sec. lIII (Jl We shall henceforth suppose that the 
width L of the lattice strip is even; the periodic boundary 
conditions in the L-direction are then compatible with 
the bipartiteness of the lattice. By virtue of the conser- 
vation law established in Sec. IIII Bl the TM has a block 
diagonal structure, with each block corresponding to a 
fixed number of particles. It is convenient to define a 
"charge" Q corresponding to each block, as Q = L/2 — n, 
where n is the number of particles. Also, we label the 
eigenvalues within each block in order of decreasing 
norm: |A?| > \\%\ > .... 

1. Correlation functions 

For entropic reasons, the largest eigenvalue must be 
located in the Q = block, A max = A". By the Perron- 
Frobcnius theorem, it corresponds to the unique eigen- 
vector in which all entries are non-negative. Consider 
first a dimer covering of a strip of size L x M, with free 



(resp. periodic) boundary conditions in the M (resp. L) 
direction. Only the TM eigenvalues of the Q = block 
will contribute to the corresponding partition function 
Z. Let us now modify the problem by marking Qo > 
vertices of the even sublattice in the bottom row, and 
Qo vertices of the odd sublattice in the top row. In the 
modified problem, dimers are required to cover all un- 
marked vertices and none of the marked vertices. The 
TM eigenvalues contributing to the modified partition 
function Zq are then exactly those of the Q = Qo block. 
(For Qo < 0, interchange the two sublattices and change 
the sign of Qo). 

Physically, the marked vertices can be interpreted as 
monomer defects in the surrounding dimer environment. 
The ratios Cq (M) = Zq q /Z define (unnormalized) cor- 
relation functions, measuring the correlations between 
the two groups of monomers, separated by a distance 
M. For M ^> L the correlations decay exponentially as 

C Qo (M)~(A?7A?) M 

If the system enjoys conformal invariance, this corre- 
sponds to an algebraic decay in the plane. More precisely, 
define the free energies per unit area as f® = L^ 1 log A^. 
The finite-size dependence 

A° - fi° = + o(L~ 2 ) (34) 

then defines a critical exponent Xq whose interpretation 
reads as follows: let Cm be a dimer covering of an N x N 
square with free boundary conditions (planar geometry) , 
with two small regions of Qo monomer defects, each re- 
gion corresponding to a definite sublattice as above. Sup- 
pose that each region has an extent of the order of the 
lattice spacing and is far from the boundaries. Then the 
probability that the two regions are separated by a dis- 
tance r satisfying 1 <C r -C N is proportional to r~ 2XQ o . 

The corresponding conformal field theory (CFT) is fur- 
ther characterized by its central charge c, which is related 
to the finite-size dependence of A max as follows [2^ : 

A° = foo + ^~ 2 + o(L- 2 ) , (35) 

where f^ — liniL^oo /" is the bulk free energy. 

While X\ determines the leading monomer-monomer 
correlation function, the leading dimcr-dimer correlation 
can be obtained from Eq. (|34|l by replacing f®° by f®. 

2. Numerical procedure 

The leading eigenvalue of a given block Q is obtained 
by an iterative procedure (the so-called power method 
j3|) in which the relevant TM block T Q is multiplied 
onto a vector of weights which is indexed by the basis 
states of that block. This vector can be taken initially as 
a single arbitrary basis state, which is known to belong 
to the block Q. The eigenvalue is then related to the 
asymptotic growth of the norm of the iterated vector. 
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L 


2 


4 


6 


8 


10 


12 


14 


16 


18 


dim(T,,) 


4 


16 


76 


384 


2004 


10672 


57628 


314368 


1728292 


dim(Ti) 


1 


8 48 


272 


1520 


8472 


47264 


264224 


1480608 


dim(T 2 ) 




1 


12 


96 


660 


4224 


26012 


156608 


929700 


dim(T 3 ) 






1 


16 


160 


1304 


9520 


65056 


426000 


dim(T 4 ) 








1 


20 


240 


2268 


18688 


141156 



TABLE I: Dimensions of the various blocks Tq of the transfer 
matrix, as functions of the strip width L. 



FIG. 8: A typical plaquette configuration. 



This procedure has multiple practical advantages: (i) 
only the iterated vector, and not Tq itself, needs to be 
stored in memory, (ii) using the factorization Eq. i|32|> 
one can take advantage of sparse matrix techniques, so 
that one iteration is performed in time ~ i.dim(Tq), 
(iii) the complete state space corresponding to Tq is au- 
tomatically generated in the iterative process. To store 
and access the weights in an efficient manner (i.e. in con- 
stant time), standard hashing techniques are employed. 

To obtain higher eigenvalues, Aj; with k > 2, one can 
similarly iterate a set of vectors which is kept mutually 
orthogonal at the end of each iteration [24]]. Alterna- 
tively, one can in some cases use the symmetry of the 
corresponding eigenvectors. As an example of this, note 
that the eigenvectors corresponding to A° (resp. A°) are 
even (resp. odd) upon shifting the lattice by one unit in 
the horizontal direction. 

The computational effort needed to obtain the largest 
eigenvalue can be judged from Table [I] which shows the 
size of the block Tq for various strip widths L. Note 
that these numbers increase much slower than the naive 
estimate 4 L , that one would obtain by considering the 
possible occupation numbers while ignoring the dimcr 
constraint and the value of Q. We limited the present 
study to T max = 18, although a couple of more sizes 
could have easily been obtained. 

The values of dim (Tq) can easily be obtained analyti- 
cally using generating function techniques. The result is 
that dim (Tq) for a given (even) value of L is the coeffi- 
cient in the term q® in the polynomial expansion of 

/ l+iq+q 2 + (l+q)^l+6q+q 2 \ 
\ ) 

^ l+4q+q 2 -(l+q)y/l+6q+q^ ^ 

The dimension dim(T) = Y1q=-l/2 °-i m (^Q) °f the total 
TM is then simply 

dim(T) = (1 + V2) L + (1 - V2) L . (37) 

This is also the dimension of the (unique block of the) 
TM when monomers are allowed; sec Sec. lVIII B 3l below. 



V. COLUMNAR ORDER AT LOW 
TEMPERATURE 

A. Possible crystalline orderings 

We present here MC results concerning the nature of 
the low-T phase. From the energy form, we expect at low- 
T a proliferation of plaqucttes containing parallel dimers. 
A natural expectation is to have a single low-T phase 
breaking the same symmetries as the ground-states in 
Fig. |21 we refer to such an order breaking both trans- 
lation and 7r/2 rotation symmetries as columnar order. 
On the other hand, from our knowledge of QDM, we 
know that another type of order could also be stabilized: 
plaquette order. We describe more precisely this order 
below. As to discriminate which kind(s) of phase(s) is 
(are) found at low-T in the dimer model Eq. (Q, we will 
introduce three different order parameters. 

1. Description of plaquette ordering 

The QDM on the square lattice is believed to have 
some plaquette long-ranged order in some finite region of 
parameter space at T = 0. In such a symmetry broken 
phase, one quarter of the square plaquettes are sponta- 
neously selected to host a pair of (quantum-mechanically) 
resonating dimers. The resulting state breaks translation 
invariance but is invariant under 7r/2 rotation with re- 
spect to the center of any plaquette (see Ref. |2f| for an 
illustration). In the quantum system, a plaquette phase 
has a (slightly) higher potential energy than a colum- 
nar crystal, but the stronger dimer resonances lower the 
plaquette state energy through the kinetic terms of the 
quantum Hamiltonian. Of course, in our classical model, 
kinetic terms are absent. Still, the thermal fluctuations 
of the dimer locations around each "flippable" plaquette 
allow to gain some entropy (compared to that of a colum- 
nar crystal) and lower the free energy. The competition 
between entropy and potential energy in the classical sys- 
tem is analogous to that between kinetic and potential 
terms in the QDM. As the plaquette phase is likely to be 
realized at T = in the QDM, it is a priori also a natu- 
ral candidate in the (finite temperature) phase diagram 
of the classical model. 

In a plaquette phase, two distant flippable plaquettes 
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are almost uncorrclatcd: if the first one is dimerized, 
say, horizontally, the second can be found in both states 
with equal probability (hence the 7r/2 rotation symme- 
try). This is not true for nearby - and necessarily cor- 
related - plaquettes, thus defining some finite correlation 
length. A columnar state can be viewed as a plaque- 
tte phase in which this plaquette-plaquette correlation 
length has grown to infinity so that all the plaquettes of 
the lattice simultaneously adopt the same orientation. A 
typical plaquette configuration is displayed in Fig.|Hl The 
plaquettc phase breaks translational symmetry but not 
7r/2-rotational symmetry. A possible scenario (eventu- 
ally ruled out by the numerical results, see below) could 
therefore be melting of the columnar crystal through an 
intermediate plaquette phase with partial restoration of 
the rotation symmetry. 



2. Order parameters 

Complex columnar order parameter — We first use the 
definition (proposed in Ref. l2fih of a complex columnar 
order parameter VPcoKr) at site r 

*col(r) = (-) r "[ft(r + x/2)-ft(r-x/2)] 

+i (-r*[n(r + y/2)-n(r-y/2)], (38) 

where x, y are unit vectors, and n is the dimer bond 
occupation number (i.e. n(r+x/2) is f if there's a dimer 
between site r and site r + x) . We define the associated 
columnar susceptibility as 



Xcoi = ^ «| £ * col (r)| 2 ) - (| £ *c ol (r)|) 2 ), (39) 



reA 



where the sums are taken only over the sublattice A. 
Another interesting quantity is the columnar Binder pT^ 
cumulant 



1*1 



2 m 



(40) 



This Binder cumulant saturates to 1/2 for a long-range 
ordered phase, and scales to in the thermodynamic 
limit for a phase with no long-range order, due to the 
Gaussian nature of the fluctuations of this order parame- 
ter. As was already noted in Ref. I2M this order parame- 
ter (and associated quantities) is sensitive to translation 
symmetry breaking and a non-zero expectation value de- 
tects both columnar and plaquette ordering. Looking at 
the phase of * co i can in principle discriminate between 
the two phases: however the phase turns out to be a noisy 
observable in our simulations and has no practical use. 
We will therefore use other indicators. 

Dimer rotation symmetry breaking — At sufficiently 
high temperatures, the system is symmetric under n/2 
rotations so that the average number of vertical dimcrs 
is equal to the average number of horizontal ones. This 



also holds in a plaquette phase, but is no longer true 
for a columnar state. A convenient way of monitoring 
the 7r/2-rotation symmetry |25j is the Dimer Symmetry 
Breaking: 



J D = 2/L 2 |A c (-)-A c (l)| 



(41) 



where N c (— ) (rcsp. N c (l)) is the number of horizontal 
(resp. vertical) dimers in the configuration c. Normal- 
ization is such that D = 1 in the pure columnar states of 
Fig. [3J It is also useful to define the corresponding sus- 
ceptibility xd = L 2 ((D 2 ) — (D) 2 ) and Binder cumulant 
B D = 1- (D 4 )/(3(D 2 ) 2 ). 

Plaquette order parameter - To discriminate posi- 
tively the plaquettc phase, we also use the following pla- 
quette order parameter 



P = 2/L 2 \J2(- 



Px+Py 



(42) 



where the sum is over all plaquettes of the lattice with 
coordinates p = (p x ,p y ), and v p = 1 if the plaquette 
with coordinates p contains two parallel dimcrs (v p = 
otherwise). This quantity can also be seen as a gener- 
alized energy at wave vector (ir,ir). The staggered fac- 
tor (— )P* + Py is essentially constant in a pure plaquette 
state, and makes the sum vanish in a columnar state. 
The expectation value (P) of the plaquette order pa- 
rameter is then in the columnar phase, and saturates 
to a finite value in the thermodynamic limit in a pla- 
quette phase. The associated plaquette susceptibility is 
Xp = L 2 ((P 2 ) — (P) 2 ) and plaquette Binder cumulant 
B P = l- (P 4 )/(3(P 2 ) 2 ). 



B. Numerical results 

1. Dimer rotation symmetry breaking 

The expectation value (D) is displayed versus T in 
Fig. Eland clearly saturates to its maximum values at low 
T. The curves for different system sizes start to differ at 
a temperature around T ~ 0.6 and in order to detect 
more finely the critical temperature T c , we use the cor- 
responding susceptibility Xd and Binder cumulant Bp. 

Xd shows a pronounced peak around T ~ 0.63 (see 
Fig. HUI) . signaling the onset of long-range order. Noticing 
that the temperature at which the susceptibility peaks 
slightly drifts when increasing system size, we differ an 
estimation of T c in favor of the Binder cumulant, which 
is known to allow accurate determinations of T c . 

The Binder cumulant Bd saturates in the thermody- 
namic limit to 2/3 at low T (see Fig. II If) and we observe 
a crossing of the curves for different system sizes for both 
cumulants at a unique temperature T c , signaling the en- 
trance into the low T columnar phase. The critical tem- 
perature is estimated from this curve to be T c = 0.65(1). 
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FIG. 9: (Color online) Dimer symmetry breaking order pa- 
rameter (D) versus temperature T for different system sizes. 
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FIG. 10: (Color online) Dimer symmetry breaking suscepti- 
bility xd versus temperature T for different system sizes. 



FIG. 11: (Color online) Dimer symmetry breaking Binder 
cumulant Bd versus temperature T for different system sizes. 




FIG. 12: (Color online) Plaquette order parameter (P) versus 
temperature T for different system sizes. The dashed line 
denotes T c — 0.65(1) as estimated by the Dimer symmetry 
breaking Binder cumulant. 



The results of this section also indicate that plaquette or- 
der is not present below T c , but leave open the possibility 
of a plaquette phase at higher T. 



2. Plaquette correlations 

The expectation value of the plaquette order parame- 
ter (P) shows a non-monotonous behavior as a function 
of T (see Fig. 1121 . with an order parameter peaking close 
to T c ~ 0.65 from above for all system sizes. One also im- 
mediately notes that (P) has overall small values and de- 
creases with system size. The plaquette susceptibility \p 
peaks slightly above T c (see Fig. ll3f) : we interpret this as 



plaquette correlations being the strongest just before the 
entrance into the columnar phase. Even though the \p 
values are very small values as compared to other typical 
susceptibilities (see for example Fig. HOf) . long-range pla- 
quette order could survive in the thermodynamic limit. 
This is clearly ruled out by the behaviour of the plaquette 
Binder cumulant (see Fig. I14|) which is non-monotonous 
as well: Bp starts to rise from its high-T zero value when 
decreasing temperature and suddenly drops down to zero 
at a temperature slightly above T c . This excludes long- 
range plaquette order. 

We conclude that (strong) plaquette correlations are 
present, start to develop as one decreases T, peak just 
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FIG. 13: (Color online) Plaquette susceptibility \p versus 
temperature T for different system sizes. The dashed line 
denotes T c = 0.65(1). 
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FIG. 15: (Color online) Columnar order parameter (I'J'coil) as 
a function of temperature T for different system sizes. 




FIG. 14: (Color online) Plaquette Binder cumulant Bp versus 
temperature T for different system sizes. The dashed line 
denotes T c = 0.65(1). 



above T c , but do not form a true thermodynamic phase 
as they are suddenly overtaken by columnar order. How- 
ever, as we show in the next section, these correlations 
possibly "pollute" the finite-size behavior of the complex 
columnar order parameter. 



Complex columnar order parameter 



The columnar order parameter 
r?(IEreA*coi(r)|) is displayed in Fig. 
ferent system sizes (from L = 16 to L = 



;i*coii) 

021 for dif- 
: 160). One 



clearly sees order setting in at low T, and the curves for 
different L start to separate roughly around T ~ 0.6. To 
determine more precisely the critical point, we inspect 
the behavior of the columnar susceptibility (Eq. (|39[) and 
Fig.HU and Binder cumulant (Eq. (gOI) and Fig. fTTI). 

The columnar susceptibility has an unexpected behav- 
ior. It exhibits two peaks: the first is quite sharp and lo- 
calized around T ~ 0.63, and the second is much broader 
around T ~ 1. It is also to be noted that whereas for 
small system sizes (L < 96), the first peak is smaller 
than the second one, the tendency is inverted for the two 
largest system sizes. This could mean that the second 
peak actually saturates to a finite value in the thermo- 
dynamic limit. Another possible scenario is that the two 
peaks merge, which is not unlikely noticing that the po- 
sitions of the maximum of the second peaks shift toward 
lower T when increasing L. Unfortunately, the currently 
available system sizes do not allow us to draw definitive 
conclusions on the scaling behaviour of Xcoi- 

The columnar Binder cumulant (Fig. I17|) also displays 
unusual features: a very flat pseudo-crossing of the curves 
for different L is observed around T ~ 1.8 (see zoom 
on left inset of Fig. I17f) and a marked anomaly around 
T ~ 0.63 (see zoom on right inset). A crossing of curves 
corresponding to different system sizes for a Binder cu- 
mulant usually denotes a transition to a long-range or- 
dered phase. However, the crossing observed here is very 
flat and it is actually almost impossible within the sta- 
tistical accuracy to locate a single crossing point for the 
three largest samples. The anomaly around T ~ 0.63 is 
particularly singular and we are not aware of such a be- 
havior being reported for a Binder cumulant in the litera- 
ture. The fact that two singularities are observed both in 
Xcoi an( l B co \ could be interpreted at first glance as signs 
of an intermediate phase. However, given our previous 
findings of strong but no long-range ordered plaquette 
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FIG. 17: (Color online) Columnar Binder cumulant B co \ as 
a function of temperature T for different system sizes. Right 
inset: Anomaly near T ~ 0.63. Left inset: Pseudo-crossing 
around T ~ 1.8. 



correlations, it appears likely that the second feature at 
high T actually disappears in the thermodynamic limit, 
whereas the first one around T ~ 0.63 subsists. Indeed, 
whereas both the columnar susceptibility and Binder cu- 
mulant are marked at the lowest temperature, the tem- 
perature at which Xcoi peaks is different from the one at 
which -B co i shows a crossing. 

We believe that the plaquette correlations found in 
Sec. IV B~2l affect the finite size behavior of other observ- 
ables and are in particular responsible for the behaviors 
observed in the columnar susceptibility and Binder cumu- 
lant. The fact that strictly speaking, xp does not peak 
where the second peak in Xcoi is present (even though the 



latter drifts with system size) might indicate that other 
types of correlations are also present above T c . 

In conclusion of this section, we find that the model 
Eq. Q shows a unique transition to a columnar order at 
T c = 0.65(1) with no intermediate phase, but with strong 
plaquettcs correlations above T c . 



VI. NATURE OF THE PHASE TRANSITION 

The model Eq.^ displaying a phase transition at T c = 
0.65(1) separating a low T columnar phase from a high T 
phase (that we will describe more carefully in Sec. IVIl| . 
we now investigate the nature of this phase transition. 



A. Energy cumulant 

We first plot in Fig.^|the behavior versus temperature 
T of the energy cumulant [2^| : 



V — 1 — 



(E 4 



3(£ 2 



(43) 



In both a disordered and in an ordered phase, this cu- 
mulant saturates to 2/3 in the thermodynamic limit |28j |. 
This is also the case at the critical point for a second 
order phase transition (even though the energy distribu- 
tion is not Gaussian), whereas for a first order transition, 
it admits a non-trivial minimum (different from 2/3) in 
the thermodynamic limit [28| . We are not aware of any 
predictions for a KT transition. 

Our results for V(T) show a clear dip close to the crit- 
ical temperature T c = 0.65(1) for all system sizes, but 
this minimum scales to 2/3 in the thermodynamic limit 
as can be clearly seen in the inset of Fig. El These results 
exclude a first order transition. 
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FIG. 19: (Color online) Specific heat per site C v /N versus 
temperature T for different system sizes. The dashed line 
indicates T c — 0.65(1). Inset: zoom on the specific heat peak. 



B. Specific Heat 

We now consider the second cumulant of the energy, 
i.e. the specific heat per site, defined as : 



a _ i (e 2 ) - {Ef 

N 



N 



J-2 



(44) 



The specific heat per site shows a pronounced peak 
which does not diverge in the thermodynamic limit (see 
Fig. ^5] and its inset) . We can see that this peak is lo- 
cated at a temperature T ~ 0.59 different from the crit- 
ical temperature T c = 0.65(1), denoted by a dashed line 
in Fig.^J For a second-order phase transition, we would 
have expected either a divergence of the specific heat at 
T c (if the critical exponent a > 0) or a cusp (for a < 0). 
Our results show that C v is completely featureless at T c : 
this is typical of a Kosterlitz-Thouless transition. We 
will see in the following sections that this is naturally ex- 
pected noticing the nature of the high-temperature phase 
that we now address. 



VII. HIGH-TEMPERATURE PHASE 

It was demonstrated more than 40 years ago that the 
T = oo point of our model (which is the classical dimer 
covering of the square lattice) is critical: the system pos- 
sesses correlations functions that decay algebraically with 
distance Q . In this section we consider the fate of various 
correlation functions in the whole high-T region ]T C , oof. 



A. Dimer-Dimer correlation functions 

We have calculated two types of dimer-dimer correla- 
tion functions in the MC simulations. Both types concern 
dimers which are chosen for simplicity in the same ori- 
entation. We take horizontal dimers without loss of gen- 
erality. The first correlator which we dub "longitudinal" 
is the connected correlation function of two horizontal 
dimers on the same row separated by a distance x : 



G\x) = (rL(r)rL(r + (z, 0))) - 1/16, 



(45) 



where n_(r) = 1 for a horizontal dimer at site r, and 
otherwise. The constant 1/16 stands for the dimer den- 
sity squared. The second one is the "transverse" correla- 
tion function of two dimers separated by a distance x on 
the same column : 



G\x) = (ra_(r)ra_(r + (0,x))) - 1/16. 



(46) 



At T = oo, the exact calculations of Ref. |j give the 
asymptotic results 



ir z x z 



(47) 



and 



G\x) 



TT 2 X 2 
1 



+ 0(x~ 3 ) 
j + 0(x- 6 ) 



x odd (48) 
x even (49) 



For all finite T > T c , we find that the longitudinal 
correlation function G l (x) remains staggered, and that it 
decays algebraically, with a decay exponent ad that varies 
continuously with the temperature T: 



G\x) ~ (-)* A{T)x- ad(T) 



(50) 



for large x (with A(T) an amplitude). In Fig. [201 wc 
represent (— ) x G l (x) for four different T = 1,2,3 and 
T = oo on a log-log scale to emphasize the power-law 
decay. The algebraic decays are eventually cut around 
L/2 due to the PBC (system size is here L = 512). The 
value of the decay exponent ad(T) can be estimated from 
these plots, however the symmetry around L/2 due to the 
PBC makes a high-precision determination of the expo- 
nent difficult, since it would depend on the range of dis- 
tances used in the fit. We will use alternative methods 
(TM calculations and winding fluctuations) in Sec. I Villi 
to estimates decay exponents. To show that the different 
ways of estimating the decay exponents are consistent, 
we have plotted in Fig. [3U| lines corresponding to the de- 
cay exponents found in Sec. lVIlil (at T = oo, we take the 
exact result otd(T = oo) = 2). These lines are in perfect 
agreement with the first part of the correlation function 
(— ) x G l (x) which is not affected by the periodicity. 

In Fig. 1211 wc plot the transverse correlation function 
G t {x) for the same T. We also find here a power-law de- 
cay with the same exponent ctd(T) as for the longitudinal 
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FIG. 20: (Color online) Staggered longitudinal dimer-dimer 
correlation function ( — ) x G l (x) versus distance x between 
dimers for four different temperatures T (log-log scale). The 
straight lines correspond to decay exponents ay(T) calculated 
in Sec. rVTTTI 



FIG. 21: (Color online) Transverse dimer-dimer correlation 
function G (x) versus distance x between dimers for four dif- 
ferent temperatures T (log- log scale). The straight lines corre- 
spond to decay exponents otd(T) calculated in Sec. IVIlTI For 
the T = oo curve, the results for even x have been omitted 
(as they are negative - see Eg. I49H . 



correlation function ((— ) X G (x) and G t {x) essentially co- 
incide for large x). Small deviations can however been 
found at small distances (see x < 10 in Fig. I21|) . where 
the data for odd or even x do not exactly coincide. This 
odd/even distinction is already present in the T = oo 
case - see Eq. (gSJ> and @SJ. We find that G t (x) is well 
fitted by the expression: 

G t (x) ~ A(T)x~ a " iT) x odd (51) 

~ A(T) X - ad{T) + B{T)x- w ^ x even (52) 

where A(T) is the same amplitude as the one found 
for the longitudinal correlation function, B(T) a negative 
constant and ui(T) a subleading correction exponent for 
even x. Our data are compatible with cu(T) ~ 2. 

B. Monomer-Monomer correlation function 



notations from M to x and Cq q to M for consistency 
with previous works). For a bipartite lattice such as the 
square lattice, M(x) = if monomers are located on 
the same sublatticc. From now on, we will consider only 
monomers on opposite sublattices. To simplify calcula- 
tions, we focus on the case x = (x, 0) (two monomers 
on the same row) . The proportionality constant is taken 
such that Af(l) = 1 HEU- 
At T = oo, the exact result M(x) ~ x^ 1 ^ 2 holds in the 
thermodynamic limit At finite T, we can estimate 
M(x) with high-precision thanks to the worm algorithm 
(see Sec. IIV A|) . Results for M{x) are displayed versus 
x on a log- log scale in Fig. [22] for the four temperatures 
used for the dimer-dimer correlation functions: we also 
observe here power-laws, with an exponent a m (T) that 
appears to vary continuously with T: 



Besides dimer-dimer correlation functions, it is use- 
ful to consider the correlator between monomers, which 
arc defined as sites not paired to any neighboring site 
by a dimer. Monomers arc by definition absent in the 
model, however we can calculate the properties of two 
test monomers in a background of dimers by considering 
the monomer-monomer correlation function 0, 0, [2|j 

M(x) - Z(x), (53) 

where Z{n) denotes the number of possible configura- 
tions where the two test monomers are separated by a 
vector x. Note that this is exactly the correlation func- 
tion C Qo (M) defined in the TM section UVBl ( we changed 



M(x) ~ x~ a ™ {T) . (54) 

For the same technical reasons as for the dimer-dimer cor- 
relation functions, we do not estimate directly from this 
plot the values of a m (T), as the error bars would be too 
large. We will obtain precise estimates of a m (T) in the 
next section lYlIII which we already display as straight 
lines in Fig.[22]to show the good agreement with the real- 
space measures of M(x). As will be understood, we note 
that whereas ctd(T) decreases when lowering T, a m (T) 
increases (this can be clearly seen with the identical color 
chart of figures I2T1 and 123)1 . 



15 




X 



FIG. 22: (Color online) Monomer-monomer correlation func- 
tion M(x) versus distance x between monomers for four differ- 
ent temperatures (T log- log scale). The straight lines corre- 
spond to decay exponents a m (T) calculated in Sec lVIlH For 
the T = oo curve, we used the exact result a m {T = oo) — 1/2. 



VIII. THEORETICAL INTERPRETATION 

We recapitulate the MC results on the finite-T behav- 
ior of the interacting dimer model Eq. with v < 0: 
we observe a transition from a high-T critical phase with 
continuously varying critical exponents to a low-T crys- 
talline phase of dimers. Even though other types of corre- 
lations are present (such as plaquette correlations), there 
is no intermediate phase in-between. The transition is 
found to be of KT type. The existence of the high-T crit- 
ical phase, the floating exponents, the KT transition and 
the dimensionality d = 2 of the problem naturally sug- 
gest Coulomb gas physics |llj |. In the following, we will 
indeed describe the (non-exact) mapping of the interact- 
ing dimer model to a CG that will account for all the 
above findings. Moreover, thanks to results of CFT and 
use of TM calculations, we will be able to determine with 
high-precision the critical exponents ctd{T) and a m (T), 
and the associated CG constant. 

The following section will also have the advantage 
of reconciliating different points of view on the model 
Eq. iJTJ. Readers familiar with height models will find 
an alternative derivation of the effective height action 
for dimer configurations. For readers familiar with CG 
physics, the TM calculations presented here will allow a 
high precision test of the CG predictions, with a precision 
probably not reachable in other models. For readers in- 
terested in the dimer enumeration problem at T — oo, an 
interesting relationship between dimer winding numbers 
fluctuations and CG constant will be derived in passing. 
Finally, for readers coming from the quantum condensed 
matter community, these results have profound implica- 
tions for the high-T regime of QDM, as will be discussed 



in Sec.lTXl 



A. Mapping to a Coulomb gas 

1. Height model 

To obtain a CG picture of the interacting dimer model, 
we first use a height description of dimer configura- 
tions ,(f. Each plaquette of the square lattice is as- 
signed a real-valued height z in the following way: going 
counterclockwise around an even (resp. odd) site on the 
square lattice, z changes by +3/4 if the bond crossed is 
occupied by a dimer and by —1/4 if it is empty (resp. 
—3/4 and +1/4). These units arc chosen such that a 
monomer on the even (resp. odd) sublattice corresponds 
to a dislocation of 1 (resp. -1) in the height. To fix the 
absolute height, one fixes the plaquette at the origin to 
have, say, z(0) = 0. By integrating out the short distance 
fluctuations of z(r), one obtains an effective action S^ff [ft] 
for the coarse-grained height h{r), defined in the contin- 
uum, which corresponds to the long- wavelength modes of 
z(r). 

Locality — The form of this effective action is con- 
strained by the fact that the microscopic model is local 
in terms of the dimer degrees of freedom. Consider a fi- 
nite area A of the system and some fixed coarse-grained 
height ft(r) for r G A. The associated free energy (ob- 
tained by summing over all microscopic dimer configura- 
tions compatible with the given h(r)) should not depend 
on the dimer positions outside A. However, by shifting 
the dimers along a closed loop, the dimer configuration 
inside is unchanged but the microscopic height z(r) is 
uniformly shifted by +1 (or —1) for all the plaquettes 
located inside the loop (and so for the coarse grained 
height h). Doing so for a large loop surrounding A, one 
therefore shows that S e s must satisfy S e g [h] = S e g [h + 1] 
for any physical h. 

7r/2 rotation symmetry — Consider a large but finite 
square area A of the lattice with linear even size L. Out- 
side A, the dimers are assumed not to cross the boundary 
of A. Let z(r) = zq + d(r) be the height inside A and zq 
the height of the, say, lower left corner of A |3(|. What- 
ever the dimer locations inside A (compatible with the 
constraint above), one can move the dimers inside A by 
a 7r/2 rotation 1Z with respect to its central plaquette. 
z is unchanged outside A but for r <E A the height is 
changed to z'(r) = zq — d(7Z(r)). Now we take a smooth 
height profile h(r € A) and evaluate the associated free 
energy S s [h] by summing over all microscopic dimcriza- 
tions giving the same coarse grained height h. By the 
rotation described above, we know that another height 
h'(r) = —h(lZ(r)) corresponds to the same set of dimer- 
izations of A, up to a rotation. Because 1Z is a sym- 
metry of the lattice and of the dimer-dimer interactions, 
both h and h' must have the same free energy and we 
get S' e ff[ft.] = S c g[—h}. Here we ignored boundary effects 
at the edge of A, which should be negligible for a large 
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FIG. 23: Assuming that the column of sites located imme- 
diately to the right side of the area A is occupied by vertical 
dimers (shaded) and that no dimer crosses the boundary of 
A (dashed line), one can shift the whole area by one lattice 
spacing to the right. The vertical dimers are then moved to 
the left side. Such a "local" translation changes the height 
profile inside A from z(r) = 20 + d(r) to z'(r) = Zo — \ — d(r). 
The change d(r) — > —d(r) reflects the fact that the two sub- 
lattices are exchanged inside A. zo is the height at the lower 
left corner of A before the translation. 



enough area. 

Translation symmetry — In addition to the assumption 
that no dimer crosses the boundary of A, we assume that 
some dimers occupy all the vertical bonds on the right 
side of A (shaded dimers in the left panel of Fig. I23fl . 
One can shift the dimers of A by one lattice spacing to 
the right provided that the column of dimers on the right 
side of A is put back on the left side after the translation 
(shaded dimers in the right panel of Fig. I23fl . If the 
height inside A is z(r) = zq + d(r), the new height after 
the translation is z'(r) — zq — | — d(r — 1). As before, z$ 
is the height of the lower left corner of A before the move 
(and Zq — j after the translation). Again, we evaluate the 
free energy S e g[h] associated to a smooth height profile 
by summing over all microscopic dimerizations giving the 
coarse-grained height h. From the "local translation" 
above, one shows that the smooth height h'(r) = — \ — 
h(r — 1) ~ —j — h{r) corresponds to the same set of 
dimerizations of A, but shifted by one lattice spacing. 
Because of the translation invariance of the model, both 
coarse-grained height profiles therefore have the same free 
energy: S c s[h] = S c s[—h — j]. As before we ignored 
boundary effects at the edge of A. We also neglected 
the variations of h at the lattice spacing scale: h(r) ~ 
h(r + 1). This is natural for a smooth height, obtained 
from a microscopic (and discrete) z(r) by filtering out 
short wavelength modes. 

From the discussion above, it appears that lattice sym- 
metries of the dimer model implies not only spatial sym- 
metries for the effective height model but also some pe- 
riodicity of the free energy in the height space. A similar 
result was previously obtained by Henley and coworkers 
using the concept of "ideal states" |3lj . 

Combining the constraints of translation and rotation 
symmetries, we get that h — > —h and h — > h + 4 should 



keep the effective action invariant. This shows that the 
only allowed "potential" terms are cos(2irph) with p an 
integer multiple of 4. As only long-wavelength modes 
(k <C 1) are kept in a coarse-grained height, only terms 
with a minimal number of space derivatives are impor- 
tant. This naturally leads to an "elastic" term (Vft.) 2 . 
We eventually get a sine- Gordon model for the coarse- 
grained height: 



dr 



7T£i|V/i(r)| 2 + ^2 V p cos(2irph(r)) 



(55) 

So far we only invoked symmetry arguments to con- 
strain the form of the effective action. It turns out that 
the elastic and potential terms with p = 4 have a sim- 
ple physical interpretation in terms of the dimer model. 
The gradient term favors "flat" heights. There arc indeed 
more dimer configurations corresponding to a flat aver- 
age height than a tilted one. This comes from the fact 
that a dimer shift along a closed loop is possible only 
if the height is constant (up to small discretization ef- 
fect) along the loop. Thus, there is more room for dimer 
moves if the average slope is small (in which case there 
are many small "iso-/i" loops available) than if the slope 
is large (fewer "iso-/i" curves). As for the cos(87r/i) terms, 
it has 4 minima (h = g , | , | , |) which precisely coincide 
with the average height of the 4 columnar configurations 
(ground-states) which minimize the dimer-dimer interac- 
tion energy (see Fig[3J). The model therefore describes a 
combination of entropy and energy (locking) effects. 

It is a standard result |13 that in 2D the relevance in 
the Renormalization Group (RG) sense of the cosine term 
depends on the period and on the stiffness constant g [33| • 
The cosine term is relevant (and locks the height) if g > 
p 2 /4 whereas it renormalizes to zero when g < p 2 /^. 
In the latter case the long-distance theory is a free field 
(elastic term only) and the system is critical ( "rough" in 
the height model terminology). The transition between 
the two phases is of the KT type. For completeness, this 
classic RG calculation is reproduced in Appendix lAl 



2. Coulomb gas 

It is well-known that the sine-Gordon model is equiv- 
alent to a (low-density) one-component CG [ll|. This 
standard mapping is reproduced in Appendix ]B\ This 
point of view has the advantage (over the sine-Gordon 
model described above) that it provides a framework to 
understand the role of the (here suppressed) monomers in 
the model, as well as other operators or correlations. 

In the mapping from the sine-Gordon model to a CG 
model, the height field is conjugate to the electric charge 
density. Similarly, it can be shown that magnetic charges 
correspond to topological defects (dislocations) in the 
height. A dual magnetic charge m = 1 (resp. m = — 1) 
corresponds to a dislocation in the height field, which can 
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be inserted "by hand" through the inclusion of monomers 
on the even (resp. odd) sublattice or through appropri- 
ate boundary conditions (see the discussion in Scc. lIVB]l . 
The exponent [34| associated to the insertion of a particle 
with electromagnetic charge (e,m) is given by |ll| 

a(e, m) = e 2 /g + gm 2 , (56) 

where g is the Coulomb gas coupling constant. The nor- 
malization of g and of the height field in Eq. (|55|) have 
been chosen such that e and m are integers and that 
standard expressions [TlT | for g are recovered. 

In the CG picture, inserting an electric charge e is im- 
plemented by a vertex operator V e {r) =: exp (2iireh(r)) : 
appearing in the Fourier expansion of any operator peri- 
odic in the coarse-grained field. Such a term with e = 1 
actually appears as a continuum limit contribution in 
the definition of the dimer operator |35j : this allows to 
identify the continuum limit of the dimer number as an 
operator with electric charge e = 1 . To magnetic charges 
(monomers) correspond dual operators h(r); however the 
fugacity for these magnetic charges is fixed to zero as we 
consider close-packed dimers. The dimer-dimer corre- 
lation function decays with an exponent related to the 
dimension of the e = 1 operator, ay — ce(l, 0) = 1/g, and 
the monomer-monomer correlation function decays with 
the exponent a m = a(0, 1) = g. This leads in particular 
to the prediction ay = l/a m , independently of T fl2T ]. 

The effective coupling constant g(T) varies with T to 
account for the continuously varying exponents. It can be 
calculated via the above mentioned relations with decay 
exponents or via fluctuations of winding numbers (see 
Scc. lVIli"H)l . As usual in a CG description, it is useful to 
have an external exact result to fix the coupling constant: 
here the calculations of Refs. 0,0 give g(T = oo) = 1/2. 

Another insight of the CG picture is about the rele- 
vance of the locking potential V p in Eq. I|55|l : the identi- 
fication with the vertex operator of an electric charge p 
immediately indicates ^l| that it will become relevant for 

2 

9 > 9c(j>) = if (this is the same result as obtained within 
the RG of Appendix El . In our model, p = 4 and g thus 
increases from 1/2 at T = oo to g(T c ) = g c (p = 4) = 4 at 
the critical point. This also means that locking potential 
terms with higher values of p will always be less relevant 
and can be ignored in the present model. 

We finally note that an equivalent route leading to the 
same CG model is to consider the low-T phase. The 
system has a four-fold ground-state degeneracy. At low 
T, a finite system is mainly composed of large domains 
where all the dimers are aligned in one of the four pos- 
sible columnar states of Fig. [3J One can associate to 
each domain orientation a q = 4 "clock spin" . One could 
therefore naively expect a transition in the 2D q = 4 state 
clock model universality class. The g-state clock model 
is well known to map to a Coulomb gas with (uncon- 
strained) integer magnetic charges and electric charges 
multiples of q [llj. The key point here is that due to 
the absence of monomers in our model, there cannot be 
isolated points (sites) around which this clock-spin can 
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FIG. 24: (Color online) Three-point fits c(L - 4, L - 2, L) 
for the central charge, as a function of W . The dashed line 
indicates T c as found by the Monte Carlo simulations. 

make a 27r rotation passing by the four possible ground 
states. This is nicely illustrated in another context in 
Ref. l36l In other words, we are considering a q = 4 
state clock model with no topological defects. Magnetic 
charges are therefore absent in the associated CG and we 
arrive at the same CG picture described above: electric 
charges can only be multiples of 4 and magnetic charges 
are absent. 



B. Transfer Matrix calculations 

TM calculations are perfectly suited to validate this 
CG scenario, thanks to the conformal invariance results 
Eq. O and PS). 

1. Fundamental exponents 

In Fig. [2] we show the finite-size estimates of the cen- 
tral charge c as a function of W = %~ V I T . The results 
shown arc three-point fits based on Eq. i|35|) with a l/I 4 
correction added. Indeed, such a non-universal correc- 
tion is predicted by conformal invariance, and is known 
to greatly improve the rate of convergence. The data 
show a clear c = 1 plateau for 1 < W < W c , where 
we estimate W c = 4.5(1) from the intersections of the 
curves. For W > W c the curves level off to zero, with the 
drop getting sharper with increasing strip width L: this 
signals the transition to non-critical behavior. The esti- 
mate of W c is in perfect agreement with the MC estimate 
of T c . 

In Fig. we show the finite-size estimates of the 
monomer exponent X\ = a m /2 as a function of W . The 
results shown are two-point fits based on Eq. (|34f> for 
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FIG. 25: (Color online) Two-point fits Xx(L - 2, L) for the 
monomer exponent, as a function of W . 




FIG. 26: (Color online) Two-point fits X 2 (L - 2, L) for the 
dimer exponent, as a function of W . 



Qo = 1 with, once again, a 1/i 4 correction added. The 
agreement with the result X\ = 1/4 (see Section Till F() 
for W = 1 is very clear. Within the critical region 
W G [1, W c ], the data can readily be extrapolated to the 
thermodynamic limit L — > oo, and one finds that X\ is a 
monotonically increasing function of W that eventually 
takes the value Xi(W c ) = 2. So the monomer perturba- 
tion is marginal at W c and relevant in the critical region, 
as predicted by the CG. For W > W c the extrapolation 
in L of the numerical data does not work well, as could 
be expected from the non-critical behavior predicted in 
this region. 

Fig. 1261 shows the finite-size estimates of the dimer ex- 
ponent X2 = ay/2 as a function of W, based on the first 



FIG. 27: (Color online) Two-point fits X(14, 16) for the first 
15 exponents in the Q = block, as functions of W . All ob- 
served exponents are simply related to the fundamental dimer 
exponent X 2 . 

two eigenvalues in the Q = block. There is an excellent 
agreement with the result X2 = 1 (see Section fill Fl) for 
W = 1. Inside the critical region W G [1, Wc], the data 
decrease monotonically, and one has to a very good pre- 
cision XiX 2 = 1/4 independently of W, confirming the 
CG scenario involving only a single coupling constant g. 
An estimate for the location of the critical point W c can 
be obtained from the crossings of the curves, and is con- 
sistent with the one given above. 

2. Higher dimer exponents 

To analyze the agreement with the CG scenario in 
more detail, we show in Fig. |23 the exponents X ob- 
tained from the first 15 gaps in the Q — block. The 
first of these coincides with the fundamental dimer expo- 
nent X2 discussed above. This level exhibits a two-fold 
degeneracy which is almost exact at finite size. The iden- 
tity operator is observed to have descendants at level 1 
and 2 (i.e. operators with constant X = 1 and X = 2), 
in agreement with CFT. 

More interestingly, all other exponents appear to be 
of the form p 2 X 2 + q, for integers p > 2 and q > 0. 
In the CG scenario, this is accounted for by operators 
of electric charge p times that of the fundamental one, 
and their descendants. Inevitably, these higher operators 
are slightly less well determined in finite size, and the 
splittings of the degeneracies which would be exact in 
the thermodynamic limit are somewhat larger. In the 
figure, we have shown in solid line style the finite-size 
data for X 2 (corresponding to p = 1), X 2 + 1 (its first 
descendant) and AX 2 (the p = 2 electric exponent). The 
agreement with the numerical data is very fine. 
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To summarize, the data of Fig. [23 provide strong ev- 
idence that the CG description of the model is correct 
and complete. 



3. Including monomers 

The generalization of the dimcr model to finite 
monomer fugacity was first studied by TM calculations 
in the grand canonical ensemble in Rcf. We give here 
a more detailed analysis. We also mention a very recent 
work by Papanikolaou, Luijten and Fradkin |l7j who in- 
dependently consider the finite doping transition in this 
model. 

The monomer exponent X\ is RG relevant {i.e. X\ < 
2) for 1 < W < W c . This means that the critical phase 
of the close-packed dimcr model is unstable toward dop- 
ing with monomers. We have verified numerically that 
there is no critical behavior for W £ [1, W c ) and finite 
monomer fugacity £. This means that the RG flow will 
be toward the trivial non-critical fixed point at £ = oo. 

However, X\ is marginal {X\ = 2) exactly at W = W c , 
and thus one may suspect, as announced in Ref. 0, that 
allowing for a finite monomer density will produce an- 
other critical line emerging from W c . To test this sus- 
picion, we have adapted the TM to accommodate for 
monomers with Boltzmann weight £. The weight of a 
pair of aligned dimers is still W = e~ v / T with v = — 1. 

Despite of this modification the basis states can still be 
described in terms of the edge occupation numbers intro- 
duced in Scc. lmGl The main difference is in the dimcr 
constraint: if a vertex is occupied by a monomer, none 
of its four incident edges may be occupied by a dimer. It 
should be noted that allowing for monomers will couple 
the different blocks Tq of the TM. Accordingly, the di- 
mension dim(T) of the modified TM is larger than in the 
pure dimer case for which one could diagonalize sector 
by sector. We have therefore restricted the study of the 
monomer-doped model to widths L < 16. The analytic 
expression for dim(T) has been given in Eq. I|37|) above. 

In Fig. EElwe show rough scans of the effective central 
charge c c ff as a function of T, for various values of £ and 
rather small sizes L. We observe that for each value of 
£, there exists a temperature T(£) for which c c ff goes 
through a maximum. The finite-size corrections to T(£) 
are found to be very small. We interpret the curve T(£) 
as a line of fixed points which correspond to the phase 
transition between high and low monomer density. 

As the maximum in c c ff becomes sharper with increas- 
ing system size (not shown), the fixed points are 
RG unstable to small variations in the parameters (£, T). 
This means that points on the low-£ side of the transition 
will renormalize toward vanishing monomer density, but 
this phase will be non-critical (crystalline phase) since 
the temperature is lower than the critical temperature in 
the pure dimer model. Points on the high-£ side of the 
transition will renormalize toward infinite £ as before. 
Note that when £ increases, T(£) decreases. This was to 




FIG. 28: (Color online) Temperature scans of the effective 
central charge c e ff(4, 6,8) for various values of the monomer 
weight £. 
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TABLE II: Phase transition temperature T(£) in the model 
where monomers are allowed with a weight £. The error bar 
on the values of T(£) is of the order ±0.0002. 



be expected, since when the dimers are diluted by more 
and more monomers, it becomes harder for them to align 
at a given temperature. 

We have determined T(£) by taking these temperature 
scans of Fig. [2^1 to larger sizes (up to L = 16) and care- 
fully studying the finite-size effects. Our final results for 
the phase transition temperatures are given in Table ITT1 

To determine whether the phase transition curve T(£) 
corresponds to a critical behavior, we show in Fig.l29lthe 
effective central charge as a function of T(£) for several 
different system sizes. We observe that there exists a 
temperature T+ such that c ~ 1 for T(£) > T*. Deter- 
mining T* from this figure is a little delicate, since even 
for the lowest value of T shown, the finite-size effects 
are such that c s decreases with increasing L. We can 
however give a first rough estimate: 



T* = 0.35 ±0.05. 



(57) 



We interpret the role of T+ as follows: For T > T*, 
the curve T(£) describes a line of second-order (contin- 
uous) phase transitions with c = 1, whereas for T < 
the phase transition becomes first-order (discontinuous). 
This means that (£*,T*) is a tricritical point. The tri- 
critical nature of this point is further corroborated by an 
easy domain-wall argument showing that the transition 
at T = is necessarily discontinuous ^3 ■ 
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FIG. 29: (Color online) The effective central charge as a func- 
tion of T"(£). The dashed line denotes c c g = 1. 
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FIG. 30: (Color online) The lowest four scaling dimensions 
x(L — 2, L) along the critical curve T(£). The dashed and 
dotted lines denote the values x = 1/8 and x = 1/2 respec- 
tively. 



To determine the universality class of the critical part 
of the curve T(£), we show in Fig. 1301 estimates for the 
lowest four scaling dimensions (denoted x\ < X2 < x% < 
Xi) as functions of T(£). A conspicuous feature is that 
x\ ~ 1/8 is almost independent of T(£), and that x-2 
is almost degenerate with X\. The existence of a con- 
stant exponent x = 1/8 (first noted in Ref. ll2T) and 
other interaction-dependent exponents, is characteristic 
of the two-dimensional Ashkin- Teller [33 (AT) model. 
We therefore conjecture that the critical curve T(£) is in 
the universality class of the AT model. 

To test this conjecture and make it more precise, we 



recall some facts about the isotropic, selfdual AT model, 
following section 12.9 of Baxter's book [H. The AT 
model is originally defined in terms of two Ising mod- 
els, each with coupling constant K, that interact via a 
four-spin (energy-energy) coupling A4. Defining lj\ = 
cxp(— K4) and L03 = cxp(A'4 — 2A"), one finds, after a 
long series of transformations, that the model is equiva- 
lent to a six- vertex model with parameters a = b = \pliO\ 
and c = V2 + ^3). The anisotropy parameter A of 
the corresponding XXZ spin chain reads 



A = 1 



1 



L03 
0J1 



1 



2cos 2 f^ 
V 4 



(58) 



where the second equality parametrizes only the critical 
regime |A| < 1, through the parameter y <E [0,2]. We 
note that y = corresponds to a KT transition. The 
original couplings, K and A4, are only real and finite for 
y € [0, 4/3), but from the point of view of the six- vertex 
model nothing special happens at y = 4/3. The critical 
exponents of the AT model differ subtly from those of 
the equivalent six-vertex model. There are two magnetic- 
type exponents, which in our notation read = 1/8 



(2) 
H 



8-4y 



and a temperature-like exponent Xt = 



and 
1 

We now claim that the critical line T(£) corresponds 
to the AT model with y £ [0, 3/2]. Moreover, we identify 

x\ = X2 of Fig.|2niwith Xfl in the AT model, X3 with x^ , 
and X4 with xt (the latter identification is only on the low- 
T side of the level crossing visible in the figure; for higher 
T we identify Xt with the analytic continuation of X4). 
The parameter y is an (unknown) increasing function of 
T, taking the values y = at T = and y = 3/2 at 
T = T C . 

Using Fig. 1301 we can now identify either from .T3 = 
1/8 or from x 4 = 1/2. These two determinations are 
consistent, but the latter leads to the best final estimate 
of the tricritical point: 



T* = 0.29 ±0.02. 



(59) 



more precise than the first estimates of Eq. \57\ and of 
Ref. [jlj For y — 0, the original couplings of the AT 
model satisfy K = A4, meaning that the two constituent 
Ising models couple strongly so as to give a single four- 
state Potts model. We conclude that the transition at 
is in the universality class of the critical ferromagnetic 
four-state Potts model. Note that this identification is 
consistent both with the tricritical nature of the transi- 
tion, and with the existence of a RG marginal direction 
in the Potts model [3?|- 

Meanwhile, for y = 3/2 we have = 1/8 and 

x H = 1/2. This is in nice agreement with the values 
of the first two electric-type exponents as obtained from 
the Coulomb gas of the pure dimer model at T = T c . 
Moreover, Xt = 2 becomes marginal and can be identi- 
fied with the monomer operator, which is responsible for 
the transition. 



21 



As a last check of the AT identification, consider the configurations are given by: 



point y = 1, where K A = and K = ^ ln(l + V2), so 
that the model decouples into two non-interacting critical 



Ising models. Here one has = 2.t^ j = 1/4 and Xt = 
1. One can verify from Fig. 1301 that this indeed happens 
at the same temperature, Ti s ; ng ~ 0.54. 

The analysis presented in this section is in agreement 
with the results of Ref . 0] ■ 



C. Winding number fluctuations and Coulomb gas 
constant 

In this section, we derive an analytical relationship be- 
tween the coupling constant g of the CG in the high-T 
region and the fluctuations of dimcr winding numbers (to 
be defined below). As winding number fluctuations are 
easily accessible in MC simulations, this allows for an in- 
dependent calculation of the CG coupling constant that 
can be compared with TM calculations. 

Orienting the dimers from the odd sublattice toward 
the even sublattice defines a fictitious "magnetic field" 
B on each bond Q. Because each site has exactly one 
dimer, the lattice divergence of this magnetic field is 
divB = 1 (resp. -1) on the even (resp. odd) sublattice. 
As a consequence (Gauss law), the flux of B through 
a contractible loop (with as many sites from each sub- 
lattice) is zero in any dimer covering |40j . However, the 
fluxes W x and W y through the two (non-contractible) cy- 
cles winding around the torus are two non-trivial integers 
In the following, we denote P(W x ,W y ) the probability 
to observe a dimer configuration with winding numbers 

W X , Wy. 

In the continuum limit, this probability can be ob- 
tained in the high-T regime where the dimer configura- 
tions are coarse-grained into a free field 0] ■ In the height 
representation of Sec. IVIII A ll the field h is defined by 
the integral of its derivative. The winding number across 
a cycle is given by this integral along the cycle. The 
probability P(W X , W y ) is the ratio of the partition func- 
tion restricted to fields having winding numbers equal to 
W x , W y across the two cycles of the torus: Zw x ,w , di- 
vided by the total partition function w gZ ^W^,W y ■ 

This ratio can be evaluated as follows: one separates 
the height field into a classical and a fluctuating part 
h = h c \ + 5 where h c \ is the solution of the equations 
of motion (a harmonic function) carrying the two wind- 
ing numbers W x ,W y , and S is a fluctuating field with 
no discontinuity. /i c i being a solution of the equations of 
motion, the crossed term disappears from the action and 
the free part of the action Eq. (|55|) is the sum of a classi- 
cal part and a fluctuating part which does not depend on 
the winding numbers. As a result, the partition function 
Zw x ,Wy factorizes into a classical part Zfy w and a fluc- 
tuating part Z' ', which being independent of the winding 
number, factorizes out of the probability P(W X , W y ). 

On a square torus of size L x ,L y , the classical height 



h(x,y) = x— +y-j—, 

L X Lly 



(60) 



and the probability is obtained by substituting this field 
in the expression of the Boltzmann weight: 



P(W x ,Wy) 



-^((T7) 2 +(^) 2 )i*L B 

„-^((xT) 2 + (f^) 2 )^ i . 



(61) 



The quadratic fluctuation of the winding number across 
the direction x is thus given by: 



J2n€Z e~ 7r9 ™ 2 ( i »/ Lx ) 



(62) 



In the high-T phase of the dimer model, short- distance 
properties are of course not described by a free height 
field. However, a local dimer move cannot change the 
winding number and we expect that the precise (non- 
Gaussian) nature of the local fluctuations will not affect 
winding number observables. (W£) is related to (fluctu- 
ations of) the global slope of the height profile. It can be 
expressed using the Fourier modes |fc| -C 1 of the height 
field only [3lJ and these are precisely those described by 
the Gaussian action. 

In the non interacting case (T = oo) where g = 1/2, 
Eq. (|62|l has recently been derived using the Pfaffian ex- 
pression of the partition function by C. Boutillier and 
B. de Tiliere . For other values of T, it allows to 
calculate through MC simulations the temperature de- 
pendence of the coupling constant g, which can be com- 
pared with the one obtained from dimer and monomer 
exponents calculated via the TM (see Sec. IVIII B l|l . 
The temperature dependence of g(T) obtained by the 
three methods match perfectly as can be seen in Fig. 1311 
where the three curves basically overlap. The value 
of g at T c (denoted by a dashed line in the figure) is 
g c = g(T c ) = 4.0(2), in agreement with the CG prediction 
g c = 4. The CG analysis and the free field calculation 
Eq. 1(:>2H are therefore entirely validated. 



IX. DISCUSSION 



Connection to other classical models 



With the interpretation of Sec. IV111 A 2l the interact- 
ing classical dimer model is naturally related to other 
scalar- field models displaying a CG behavior ^lj. Even 
if it is not probably exactly solvable, its definition puts 
it among the most simple ones, along with the XY or 
6-vertex models. Moreover, it can be quite naturally ex- 
tended (e.g. by changing the sign of interactions, going 
to other geometries), which allows to test various features 
of the CG. 

The model Eq. also admits a height description 
(Sec. IVIII A l|) . which relates it to a variety of spin 
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FIG. 31: (Color online) Coulomb gas coupling constant g 
as a function of temperature T obtained via TM (through 
exponents Xi and X2 - see Sec. lVIII B 1(1 and MC calculations 
(thanks to Eq.Id^. 



and Solid-On-Solid models (see Ref. 03 and references 
therein) which display the same physical behavior. Wc 
believe this interacting dimer model offers several advan- 
tages over these models: (i) its definition is quite sim- 
ple and natural, (ii) it admits two exacts points (at zero 
and infinite temperature), (iii) very efficient numerical 
simulations are possible (thanks to the directed-loop al- 
gorithm 0| for the MC simulations and to the relative 
small size of the configuration space for the TM calcula- 
tions). 

Concerning numerical simulations of height-like mod- 
els, we find that the height description is not strictly nec- 
essary to obtain reliable numerical results (as opposed to 
Ref. l3ll c) . Henley and coworkers [3lJ find the height rep- 
resentation particularly useful because it allows for exam- 
ple to calculate the stiffness (or CG) constant g via the 
long-wave lengths fluctuations of the Fourier transform 
of the height field. We have proposed here an alternative 
method to obtain g via the fluctuations of the winding 
number (see Sec. rVIliT!|) . which has the advantage of not 
involving any fitting procedure. 

Concerning the physics of classical dimers, the anal- 
ysis presented here can probably be extended to the 
recent results of Sandvik and Moessner 0| on non- 
interacting models with longer-range dimers. For models 
with dimers allowed between next-nearest neighboring 
sites, the bipartite nature of the lattice is lost and the 
(dislocation-free) height mapping no longer valid: cor- 
relations become exponential. For models with a frac- 
tion of dimers allowed between fourth-nearest-neighbors, 
both bipartiteness and height descriptions are present: 
the system stays critical. There, Ref. finds that the 
monomer-monomer decay exponent continuously varies 
(from a m = 1/2 to a m = 1/9) with the fraction (fu- 



gacity W4) of longer-range dimers: this is naturally ex- 
plained in our framework by considering the variation of 
the CG constant g with W4 - which in turns controls the 
monomer-monomer decay exponent. It it also found that 
the dimer-dimer correlation keeps its 1/r 2 behavior. This 
last finding was accounted for in Ref. by considering 
the "dipolar" terms in the dimer-dimer correlation func- 
tions. The corresponding part of the dimer-dimer cor- 
relation function continues to vary as r~ 2 whereas the 
vertex contribution [^3 scales as r -1 /^™ 4 ). The long- 
distance behavior observed in the numerics is therefore 
dominated by the largest exponent between both; here it 
is 2. This is reminiscent of spin-spin correlation functions 
of the XXZ quantum spin chain ^3] where the uniform 
part of the (S z (0)S z (r)) correlation function decays as 
1/r 2 whereas the staggered part decays as l/r-^ A ) where 
/ is a continuous function of the anisotropy parameter 
A. 

Finally, we remind that the interacting dimer model 
Eq. Q was introduced in the context of liquid crystal 
physics. By means of a low-T expansion, Poland and 
Swaminathan estimated in Ref. Il3b the transition tem- 
perature to be T c = 0.61(1), quite close to the actual re- 
sult. However, they incorrectly stated at that time that 
the transition was Ising-like (second order) and did not 
realize the critical nature of the high-T phase. 



B. Implications for finite temperature properties of 
the Quantum Dimer Model on the square lattice 

The classical model naturally offers informations about 
the finite-T properties of QDM @, where quantum ef- 
fects are not dominant. In particular, the critical phase 
found in our model should be present in all the high-T 
phase diagram of the QDM: indeed, the rough nature of 
this phase was shown to be intimately tied to the dimer 
hard-core constraint and not to the details of fluctuations 
(thermal or quantum). Wc note in passing that, as in the 
classical case, doping the QDM with monomers (static or 
mobile) will immediately destroy this critical phase [2flj . 
We also think that the finite-T melting of the columnar 
phase will proceed via a KT transition as described here 
everywhere in the phase diagram of the QDM (presum- 
ably, this is also the case for the melting of the plaquette 
phase) . We note that these findings are in contrast to the 
finite-T phase diag ram of the QDM proposed by Leung 
and coworkers |25| . who speculated that the columnar 
crystal would melt in two steps (with an intermediate 
plaquette phase), even in the classical limit. 

The T = phase diagram of the QDM is usually ac- 
cepted as this [25L I44I l4Sl | (with the standard notation 
of t for the kinetic energy gained by flipping a plaque- 
tte): for large negative v/t, the system is in a columnar 
phase. Increasing v/t, the QDM experiences a quantum 
phase transition to a plaquette phase at a quite uncer- 
tain value of v/t (see discussion below). The plaque- 
tte phase ends up at the Rokhsar-Kivelson (RK) multi- 
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v/t 



Columnar Plaquette RK Staggered 



FIG. 32: (Color online) Schematic finite-temperature phase 
diagram of the Quantum Dimer Model on the square lattice. 
The different crystals (columnar, plaquette and staggered) all 
melt into a high temperature critical phase (denoted in blue) 
with algebraic correlations characterized by a Coulomb gas 
constant g. See text for details. 



critical point v/t = 1. where the correlation functions 
display algebraic quasi-long range order |j| . For v/t > 1 , 
the system is in the staggered phase. 

We are now armed to connect this phase diagram to 
our finite- T results for t = and draw our proposed 
(v/t, T/\v\) schematic phase diagram (see Fig. I32|l . We 
believe the critical phase present in the whole high-T re- 
gion extends down to low temperatures at the RK point. 
This critical phase can be parametrized by the Coulomb 
gas coupling constant g, and we draw lines of iso-g in 
this region. The line g = 4 denotes the boundary be- 
tween the high-T phase and the columnar phase (and 
presumably to the plaquette phase as well), correspond- 
ing to the expected KT transition. As the transition 
from the plaquette to the columnar phase is expected on 
general grounds to be first order at T = 0, we expect 
the transition to subsist at finite temperature. On the 
right hand side of the phase diagram, the rough phase 
will give rise to the staggered phase at low temperatures 
via a vanishing of the CG constant g — > 0. First results 
indicate that the transition is here first order (IH |46T ]. 
Finally, we conjecture that all iso-g lines meet at the RK 
point (T = 0, v/t = 1) (see Fig. 13211 . consistent with the 
multi-critical nature of the RK point jj^j . Even though 
numerically difficult (see discussion below), it would be 
very interesting to test these predictions, especially in the 
vicinity of the RK point. 

We now make a remark on the existence of the pla- 
quette phase. We first note that our classical model is 
amenable to large-scale simulations, which is not the case 
of the QDM. Numerical simulations of the QDM are in- 
deed very difficult: exact diagonalizations |25l |26[ are 



limited to small sizes (at most 10 x 10 lattices) and even if 
the sign problem is absent in the QDM, Quantum Monte 
Carlo (QMC) calculations are notoriously hard. Some 
progresses have however recently been made in the Green 
Function QMC formulation |4J, |45j, |4jfl . These difficul- 
ties could have important consequences for the plaquette 
phase: indeed, we saw that plaquette correlations were 
important in our model and altered finite-size behavior, 
even on systems as large as 160 x 160. This is an in- 
dication that the plaquette phase in the QDM could as 
well just disappear in the thermodynamic limit. This 
is not so unlikely noticing that the critical point sepa- 
rating the columnar and plaquette phase was first re- 
ported [2f| to be at v/t ~ —0.2, and then presumably 
around v/t = 0.60(5) when larger samples were avail- 
able [44I l45l| : this indicates that the plaquette phase ex- 
tent actually shrinks when increasing system size. 



X. CONCLUSION 

In this work, we studied a model of interacting close- 
packed dimers on the square lattice, where the interac- 
tion tends to align neighboring dimers. The ground-state 
constitutes of dimers aligned in column and is four-fold 
degenerate. The T — 00 point is the dimer coverings enu- 
merating problem, and displays algebraic correlations. 
With the help of MC and TM simulations, we could show 
that the low-T crystalline columnar phase leaves place 
to a high-T critical phase, with floating decay exponents 
of the correlation functions. The transition takes place 
at T c = 0.65(1) and is of the KT type. Via a height 
description of the dimer configurations, we determined 
the effective continuum theory with a competition be- 
tween a stiffness term which favors rough height profiles 
(critical phase for the dimers), and a locking potential 
favoring heights to be locked on four particular values 
(forcing the dimers to align in one of the four colum- 
nar ground states). This description is equivalent to a 
CG of electric charges, and the locking potential is in 
this picture associated with the vertex operator of e = 4 
electric charges. This allows us to interpret all the nu- 
merical results obtained in the high-T phase in terms 
of the CG coupling constant g, which we calculate nu- 
merically with high precision. The transition is under- 
stood in the CG language as a proliferation of e = 4 elec- 
tric charges. This model is probably one of the simplest 
and versions of a CG with floating exponents. Finally, 
we have established that doping the dimer model with 
monomers produces another critical line emanating from 
the KT point. There is strong numerical evidence that 
this line is in the universality class of the Ashkin- Teller 
model. The AT line terminates in a tricritical point at 
finite T* = 0.29(2), in the four-state Potts universality 
class. Below T* the transition goes first-order. Besides 
their relevance to other classical models, the results of 
this work also have important implications for QDM on 
bipartite lattices (see Sec. lIXf) . In particular, it indicates 
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that their high-T phase is also critical, and we conjecture 
that this criticality extends down to the RK point. 
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APPENDIX A: PERTURB ATIVE 
RENORMALIZATION GROUP FOR THE 2D 
SINE-GORDON MODEL 

We follow here the presentation given by Levitov|4S 
We write the energy of the sine-Gordon model as: 



E 



i(V$) 2 + Acos(/3$) 



(Al) 



We introduce two cut-offs A and A' in momentum 
space so that the field <f> splits into fast and slow compo- 
nents: 



$(fc) = 



(A2) 



where $(fc) > (fast) is non-zero only if A' < \k\ < A and 
&(k) < (slow) is non-zero only if < |fc| < A'. 

As usual we integrate over the fast component < I> > to 
obtain a rcnormalizcd energy for the slow degrees of free- 
dom: 



e = / p[d>>] e - £ (* <+$ ^. (A3) 



It is simple to see that the elastic part of the energy does 
not couple the slow and fast components 



so that we have: 



(A4) 



D -E A (*<) 



-i/d D r(V* < 



,>\2 



In the limit where the cosine term can be treated pertur- 
batively we have: 



E (*<) + § / d"r(V* < ) 2 



X H( 1 + Acos(/3($>(r) + $<(r)))) 

r 

~ n( i+A ( c ° s (/ 3 ( $> w+ <i,< w)))*>) 

r 

~ exp^y d D rco8(f3($>(r) +$<( r )))^ 
and finally: 
E A \<S><) ~ i J d D r{^^<) 2 

+X{ [ d D rcoa(p(^> > (r) + *<(»•))) )(A7) 



where (•••)$> is an expectation value with respect to 
the (Gaussian) weight e~i/ d r ( v * > ) . This expectation 
value can be computed in the following way: 



and 



<cos(/?(<i>>(r) + <I><(r)))) $> = 

l e ^* < (r) / e H34>>(r)\ +R 

2 \ /$> 



(A8) 



<|fc|<A 



d D k 

(2tt) d 



The last expression is a product of Gaussian integrals. 
The result is: 



3 i j 3* > (r) 



d 2 



~ exp 
*> \ 2 



d D k 1 



|(A10) 



/A'<|fc|<A W) D k 2 J 

Combining this result with Eqs. (|A7|) and l|A8() we get: 
E A '($<) ~ 1/ d D r(V$<) 2 

+\* J d D rcos(/3($<(r))) (All) 

with 

(3 2 [ d D k 1 " 



A* = A exp 



If D = 2 we have: 



A'<jfc|<A 



(2tt) d fc 2 



(A12) 



x J V [$>] exp \ -\ J d D r{V& 
+A / d D rco S (/3($> + $<))) . (A5) and thus 



d 2 k 1 1 , A 



A'<| fc |<A (2tt) 2 k 2 2tt °A 



log- (A13) 



A* = A 



A , N P 2 /(4tt) 



(A14) 



The elastic term is of the order of A 2 for and of 
order A' 2 for $ < . By integrating out <!>>, the importance 
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of the elastic term has been reduced by a factor (x) ■ 
The strength of the cosine potential relative to the elastic 
terms has thus been multiplied by 



(A15) 



when going from E to E A . Thus, we find that the cosine 
term is irrelevant if /3 2 > 8ir. On the other hand, the 
systems goes into a locked phase when /3 2 < 87r. 

To make contact with Eq. I|55|). we take the following 
normalization: 



E= I d D r irg (V/i) z + V p cos(27rp h) 



(A16) 



and we find that it is equivalent to Eq. i|Al|) provided 



(3 = Py With this normalization, the cosine term 

is relevant when g > p 2 /4. For p = 4 (relevant for the 
square lattice dimcr model), we obtain that the cosine is 
relevant when the stiffness g exceeds g c = 4. 



APPENDIX B: MAPPING OF THE 
ONE-COMPONENT COULOMB GAS TO THE 
SINE-GORDON MODEL 

This derivation is a slightly more detailed version of the 
argument presented in chapter 4 of Polyakov's book |50j| . 
We start from a system of charges interacting with a 
Coulomb potential V(r): 



Z = E^7[EE cx p(-E 1aVbV(r a 

N '{«a}K} V a^b 



n) 



(Bl) 



where N is the number of charges, £ their fugacity, r a 
(a = 1 • • • N) their coordinates on a D-dimensional cu- 
bic lattice (with unit lattice spacing) and q a £ Z their 
charges. 

We write the interaction in momentum space: 

^2l-1bV(r a -r b ) = $>(r)?(r / )V(r-r / ) (B2) 

a^b r,r' 



d D k 

(2tt) d 



where 



q(r)=^q a S(r 



q(k)q(-k)V(k)(B3) 



(B4) 



is the charge density at r and the Coulomb potential is 



V(k) 



K 
ifc 2 " 



(B5) 



We decouple the charges (Hubbard-Stratonovich) by in- 
troducing a real scalar field x( r ) : 



Z = E^EE /^xM]cxp(-l]T(Vx(r)) 2 

N ' {la} {r a } \ r 



Then we will perform the summation over the particle 
degrees of freedom: their number N, their charges q a and 
locations r a . This is just a sum over the charge density 
q(r): 



E 



E 

Mr)} 



(B7) 



In the limit of small fugacity £ <C 1, configurations with 
\q(r)\ > 1 are exponentially suppressed, and we can keep 
only q(r) E {—1,0, 1}. The number of particles is thus 
JV = 2~)r ( ?( r ) 2 - We have to evaluate: 



E ^ 



l{r) p iJ2 r X(r)q(r) 



9 (r)e{-l,0,l} 



=n E^ 2gix(r)9 

r \q=-l ) 

= H(l + 2tco8 X (r)) 

r 

exp ^2^^cosx(r)^ . 



(B8) 

(B9) 
(BIO) 

(Bll) 



In the limit £ <c 1, the partition function is that of the 
sine- Gordon model: 



Z~J p[ x ( r )] CX p^--i^(Vx(r)) 2 + 2e^cos X (r)^ 

(B12) 
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